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ABSTRACT 

We study the efficiency and dynamics of supermassive black hole binary mergers driven 
by angular momentum loss to small-scale gas discs. Such binaries form after major 
galaxy mergers, but their fate is unclear since hardening through stellar scattering 
becomes very inefficient at sub-parsec distances. Gas discs may dominate binary dy- 
namics on these scales, and promote mergers. Using numerical simulations, we in- 
vestigate the evolution of the semi-major axis and eccentricity of binaries embedded 
within geometrically thin gas discs. Our simulations directly resolve angular momen- 
tum transport within the disc, which at the radii of interest is likely dominated by 
disc self-gravity. We show that the binary decays at a rate which is in good agree- 
ment with analytical estimates, while the eccentricity grows. Saturation of eccentricity 
growth is not observed up to values e ^ 0.35. Accretion onto the black holes is vari- 
able, and is roughly modulated by the binary orbital frequency. Scaling our results, 
we analytically estimate the maximum rate of binary decay that is possible without 
fragmentation occuring within the surrounding gas disc, and compare that rate to 
an estimate of the stellar dynamical hardening rate. For binary masses in the range 
10 5 M Q ;$ M 10 s M Q we find that decay due to gas discs may dominate for sepa- 
rations below a ~ 10~ 2 pc — 0.1 pc, in the regime where the disc is optically thick. 
The minimum merger time scale is shorter than the Hubble time for M ;$ 1O 7 M . 
This implies that gas discs could commonly attend relatively low mass black hole 
mergers, and that a significant population of binaries might exist at separations of a 
few hundredths of a pc, where the combined decay rate is slowest. For more massive 
binaries, where this mechanism fails to act quickly enough, we suggest that scattering 
of stars formed within a fragmenting gas disc could act as a significant additional sink 
of binary angular momentum. 

Key words: black hole physics - accretion: accretion discs - galaxies: nuclei - galax- 
ies: active 



1 INTRODUCTION 

The assembly of present-day galaxies occurs via the hierar- 
chical merger of lower mass progenitors, many of which seem 
likely to have supermassive black holes (SMBHs) in their nu- 
clei. Following a merger of two galaxies, each harbouring a 
black hole, we expect the SMBHs to sink toward the centre 
of the merger product and eventually form a black hole bi- 
nary. What happens subsequently is less clear. There is no 
doubt that at sufficiently small separations (a £ 10~ 3 pc) 
loss of angular momentum to gravitational radiation occurs 
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rapidly enough to effect coalescence. However, at larger sep- 
arations an as-yet undetermined mix of angular momentum 
loss to stars and to gas is required to bring the black holes 
into the gravitational radiation inspiral regime. Observa- 
tions provide scant information. Detecting black hole bina- 
ries on pc or sub-pc scales (where the black holes' gravity 
dominates the galactic potential) is extremely difficult, with 
the closest confirmed binary in the radio galaxy 0402+379 
(Rodriguez et al., 2006) having a projected separation of 
7.3 pc. A variety of circumstantial evidence - of which the 
most compelling is probably the small observed scatter in 
the relationships between black hole mass and galaxy prop- 
erties (Ferrarese et al., 2006; Gebhardt et al., 2000) - sug- 
gests that binaries do merge, though these arguments do 
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not constrain the time scale or the mechanism that facili- 
tates mergers. 

One longstanding motivation for considering the role 
of gas in hastening binary mergers has come from the fact 
that stellar dynamical processes may simply fail. Although 
dynamical friction against the stellar background is rapid 
for large separations (Milosavljevic & Merritt, 2003), for 
pc-scale binaries the rate slows dramatically as the black 
holes eject the finite number of stars that have orbits that 
come close enough to the centre to interact with the bi- 
nary. Replenishment of stars on these loss cone orbits occurs 
on the relaxation time, which typically exceeds the Hub- 
ble time (Begelman et al., 1980; Yu, 2002; Merritt, 2006). 
This 'last parsec' problem occurs for spherical or axisym- 
metric nuclei, and is particularly severe for massive galax- 
ies which are observed to possess rather low density stellar 
cores (Faber et al., 1997). More recently, interest in this 
long standing problem has been further piqued by the real- 
ization that if gas is responsible for driving mergers, then it 
is likely that some of that material will survive in the im- 
mediate vicinity of the holes at the moment of coalescence. 
Plausibly, some small fraction of the energy released dur- 
ing the merger may go into heating the gas, producing an 
electromagnetic precursor (Armitage & Natarajan, 2002) or 
afterglow (Milosavljevic & Phinney, 2005; Lippai et al., 2008; 
Schnittman & Krolik, 2008; Shields & Bonning, 2008) to the 
gravitational wave event. Detection of an electromagnetic 
counterpart would greatly improve the accuracy with which 
space-based gravitational wave observatories such as LISA 
can localize mergers, and allow much more information to 
be gathered about the properties of the host galaxy (Kocsis 
et al., 2007). It may also be possible to determine whether 
gas discs typically attend mergers by searching for periodic 
electromagnetic signals produced by binaries at larger sepa- 
rations, independent of any gravitational wave information 
Haiman et al. (2008). 

Theoretically, it is well established that the dissipative 
nature of gas allows inflow down to pc scales in galactic 
nuclei following galactic mergers (Mihos & Hernquist, 1996; 
Springel et al., 2005; Mayer et al., 2007; Levine et al., 2008). 
Perhaps more surprisingly, smaller amounts of gas appear 
to be able to penetrate deep within the black hole's sphere 
of influence even in galaxies such as the Milky Way that are 
far removed from any significant merger activity. Evidence 
for such inflows comes from observations of young stars with 
disc-like kinematics close to Sgr A* (Paumard et al., 2006; 
Lu et al., 2008), which may well have formed in situ from 
gas at sub-pc scales (Levin & Beloborodov, 2003; Nayakshin 
& Cuadra, 2005). Taken together, these lines of argument 
suggest that it is plausible to invoke the existence of gas 
discs at small radii as a common feature of galactic nuclei, 
though the typical masses and structure of such discs are 
subject to considerable uncertainty. Prior theoretical work 
has established that gas discs are likely to be important for 
SMBH binary mergers, provided that the mass of gas in the 
disc is roughly of the same order of magnitude as the mass 
in black holes (Gould & Rix, 2000; Armitage & Natarajan, 
2002; Escala et al., 2005; Dotti et al., 2007; MacFadyen & 
Milosavljevic, 2008; Hayasaki, 2008). 

In this paper we revisit the interaction between a SMBH 
binary and a relatively low mass circumbinary gas disc. Our 
main innovation is to study this interaction using numerical 



simulations that directly resolve the physical mechanism of 
angular momentum transport within the gas. At the radii of 
greatest interest for the SMBH binary merger problem (typ- 
ically tenths of a pc, where the stellar dynamical time scales 
are longest for many galaxies) it is expected that self-gravity 
will dominate the dynamics of surrounding gas discs (Shlos- 
man et al., 1990; Goodman, 2003). The onset of self-gravity 
in an accretion disc can result either in a quasi-stable disc 
with outward angular momentum transport, or in fragmen- 
tation of the disc into stars. Either possibility is of great 
interest for its effect of SMBH binary mergers. Here, we 
study directly the first regime in which the disc remains 
stable and angular momentum is transported by the action 
of self-gravitating spiral arms. Physically this occurs if the 
cooling time in the disc exceeds the local dynamical time 
(Gammie, 2001). We use the simulations to study the rate 
at which the gas drives the binary toward merger, and the 
effect that the gas disc has on the eccentricity of the bi- 
nary. Previous simulations that have used either artificial 
viscosity or a Navier-Stokes formulation 1 to model angular 
momentum transport have found that eccentricity growth 
(both of the binary, and within the gas disc) attends the de- 
cay of the binary semi-major axis (Artymowicz et al., 1991; 
Armitage & Natarajan, 2005; MacFadyen &: Milosavljevic, 
2008). The time dependent accretion signatures that result 
may permit the identification of black hole binaries within 
gas rich systems prior to the final coalescence phase. 



2 INITIAL CONDITIONS AND NUMERICAL 
METHOD 

In this paper we present simulations that follow the evolu- 
tion of a binary that initially has a total mass M, separation 
ao, and orbital frequency D.q = ^/GM/afy The mass ratio 
between the two components is set to q = M1/M2 = 3, 
as typically expected for a major galaxy merger (Volonteri 
et al., 2003). We try two different initial orbital configura- 
tions for the binary: a circular orbit of radius ao and an 
eccentric orbit with semi-major axis size ao and eccentricity 
e = 0.1. 

Around the binary we set a co-rotating circular disc 
of gas. The disc is aligned with the orbital angular mo- 
mentum of the binary, has a total mass of Md = 0.2M 
and extends initially from 2ao to 5ao. The surface den- 
sity of the disc decreases with radius R (measured from 
the binary centre of mass) as S(J?) oc R~ . The disc di- 
mensionless thickness is initially constant H/R = 2Md/M. 
We set the initial internal energy of the gas accordingly as 
u = (H/Rfvi/i-yi-y - 1)), where v K = \/GW/R is the Ke- 
plerian velocity around the binary and 7 = 5/3 is the adi- 

1 Although the use of a Navier-Stokes viscosity is an improve- 
ment over a reliance on purely numerical effects to transport an- 
gular momentum, it too may lead to unphysical behaviour. In 
particular, it is known that a disc subject to only a Navier-Stokes 
shear viscosity is unstable to the growth of eccentricity even in 
the absence of perturbations (Ogilvic, 2001). In general, there is 
no reason to assume that a 'turbulent viscosity' resulting from 
a physical process such as self-gravity or the magnetorotational 
instability (Balbus & Hawley, 1998; Balbus & Papaloizou, 1999) 
will behave in the same way as a microscopic fluid viscosity. 
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Figure 1. Logarithmic maps of the disc column density (in units of Ma ) at different times during the simulation. The panel at t = 
shows the smooth initial conditions. Until t sa 350f2^ , material piles up at R ~ 3ao as a result of the torques shown in Fig. 3, forming 
a dense ring. The ring breaks due to its self-gravity (see Fig. 4), spreading the gas approximately over the original radial range. A spiral 
pattern develops, and the disc stays in that state until at least t as 1200S!q , when the simulation ends. 



abatic index of the gas. The disc stability to self-gravity is 
measured by the Toomre parameter Q = Qcs/^CE), which 
translates to Q w (H/R)(M/Md) for a disc in hydrostatic 
equilibrium. In our case, the thickness of the disc ensures 
that the disc is not initially unstable to self-gravity, as the 
Toomre parameter is « 2. Only after cooling has affected 
the gas, will the disc get thinner and become unstable. 

To calculate the evolution of the system we use the 
smoothed particle hydrodynamics (SPH; e.g., Monaghan, 
1992) code Gadget-2 (Springel, 2005). The code solves for 
the adiabatic hydrodynamics and gravitational forces of the 



system, and includes a term for artificial viscosity, necessary 
to treat shocks. 

On top of the basic hydrodynamics set-up we add cool- 
ing. This is implemented by defining a cooling time and 
setting the radiative cooling term for each gas particle to be 
(du/dt) cooling = — w/tcooi. The cooling time is set to be pro- 
portional to the orbital time of the gas around the central 
binary, t coo \ = /3/fi, with j3 a constant and O. = \J GM/R 3 . 
This prescription is commonly used in the literature (e.g., 
Gammie, 2001; Rice et al., 2005; Nayakshin et al., 2007; 
Alexander et al., 2008), since the condition for disc frag- 
mentation corresponds to a constant /3. As we will discuss 
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later, in a real disc changes to the surface density and opac- 
ity with radius result in a value of f3 that varies with radius 
- we include this physics in our analytic estimates but not 
in the numerical work. In this paper we concentrate on the 
case where the cooling is not fast enough to make the disc 
fragment, therefore we set the value of j3 = 10. We defer 
the regime of faster cooling and subsequent star formation, 
P ^ 5, to future study. 

The gravity among gas particles is calculated using a 
Barnes-Hut tree, as implemented in Gadget-2. However, 
since higher accuracy is required when calculating the evo- 
lution of the binary orbit we compute the gravitational 
forces on the SMBHs directly, by summing exactly the forces 
from all gas particles instead of using a tree approximation. 
Maintaining symmetry, the gravitational attraction from the 
SMBHs is added directly to each particle. 

While most of the gas is expected to remain in the disc, 
accretion streams allow matter to flow from the inner edge 
of the disc onto the black holes, where it becomes tightly 
bound. To prevent the computation from halting due to the 
short time-steps required to follow these particles, we 'ac- 
crete' all gas that gets within a short distance (a sink radius 
of O.lao) of either black hole (Cuadra et al., 2006). The mass 
and momentum of the accreted gas are added to the respec- 
tive black hole, thereby conserving the linear and angular 
momentum of the system 2 . When we discuss the 'accretion 
rate' onto the black holes, the reader should be aware that it 
is this accretion of particles that pass within the sink radius 
that is the actual quantity being computed. The physical 
accretion rate would track the numerical one provided that 
the residence time of gas through the (unmodeled) disc at 
smaller radii is shorter than the characteristic time scale of 
accretion variations. 

We run each calculation for more than a thousand bi- 
nary dynamical times using 2 million SPH particles to model 
the gas. In addition, to test the numerical convergence of our 
results, we run shorter simulations with 8 million particles. 



3 EVOLUTION OF THE DISC 

All simulations display very similar disc evolution, so in this 
section we concentrate on the run with 2 million gas particles 
and the binary on an initially circular orbit. Figures 1 and 
2 show the evolution of the surface density E of the disc. 
Figure 1 shows maps of E at 6 selected times, while Fig. 2 
shows the azimuthally-averaged E(i?) profiles at the initial 
conditions and every lOOf^ 1 in the interval 50 < if2o < 550. 

The disc has initially a density profile E(i?) oc R~ l , as 
seen in the solid line of Fig. 2. The sharp boundaries are 
quickly relaxed and become smoother. Figure 3 shows the 
different contributions to the torque per unit radius acting 
at this early stage in the simulation. The gravitational po- 
tential of the binary produces a torque that oscillates as a 
function of radius, as described in detail by MacFadyen & 
Milosavljevic (2008). On top of that there is a large compo- 
nent produced by the disc self-gravity, and smaller contri- 
butions from gas pressure and artificial viscosity. The total 



2 This has been verified a posteriori. The change in total mo- 
mentum of the system is only fewx 10 — 4 Mao£2o after lOOOSl^ 1 , 
accurate enough for our purposes. 
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Figure 2. Azimuthally-averaged column density of the disc at 
times tn = 0,50, 150, 250... 550. Initially (solid line), the disc 
has a surface density profile £ oc R~ 1 and extends for 2 < 
R/ao < 5. Until t ^ 350Cq (dotted lines) the material piles up 
at R Kb 3ao and forms an increasingly dense ring. This ring finally 
becomes unstable due to its self-gravity and breaks, making the 
density profile flatter for t ^ SSOf^ 1 (dashed lines). 




Figure 3. Azimuthally-averaged differential torques dG/dR av- 
eraged over the interval ff!o =30-50. Torque is dominated by 
gravitational effects, both due to the binary (green line) and the 
disc itself (purple). Pressure (blue) and artificial viscosity (red) 
effects are negligible, except at the disc edges, where a very small 
amount of gas is affected. 



torque has a minimum at R/ao s=s 3, where material piles 
up and eventually a dense ring forms. There are other min- 
ima in dG/dR(R) - notably one at R/ao w 2 - but they 
occur at locations where not much gas is present, prevent- 
ing the formation of other density enhancements. From the 
figure it is clear that the torques produced by hydrodynam- 
ical forces (both pressure forces, which could be physically 
important, and artificial viscosity which is a numerical con- 
cern) are negligible compared with the gravitational torques. 
The artificial viscosity torque becomes somewhat important 
only at the edges of the disc, due to the fact that shocks 
form at these locations. 

As seen in Fig. 2, the gas piles up at R « 3ao, until the 
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Figure 4. Azimuthally-averaged Toomre parameter of the disc at 
different times. The different lines correspond to the same times 
shown in Fig. 2. The disc is initially stable to self-gravity, as Q > 1 
(solid line), but this value at R £3 3ao decreases with time as gas 
accumulates in a ring at that location (dotted lines). After the 
ring breaks, Q remains larger than unity (dashed lines). 



column density becomes so high that the Toomre parame- 
ter Q = S2c a /(7rG£), plotted azimuthally averaged in Fig. 4, 
becomes smaller than unity and the ring collapses due to 
its self-gravity. That process can be see in the panels 2-4 of 
Fig. 1, where the ring forms, develops asymmetric structure, 
and finally collapses after reaching values Q ~ 0.5. The gas 
that was in the ring then disperses in radius and reaches 
a quasi-steady state characterised by a complicated spiral 
pattern and a value of Q m 1.5. It is interesting to note that 
the critical value of Q is different when the gas accumulates 
in a ring, at the early stage of the simulations, compared to 
the case where the gas is spread over a larger radial range, 
later in the evolution of the disc. This shows that the disc 
behaviour is not completely determined by its (azimuthally 
averaged) local properties. Instead, its stability seems to de- 
pend on the properties of the gas at neighbouring annuli and 
on its azimuthal distribution (see, e.g., Lodato, 2007). 

To characterise the spiral structure that forms in the 
disc, we measure its Fourier modes, defined as 



RdR\ d0E(7?,0)exp(-im0)|. (1) 
Jo 



Figure 5 shows the evolution of the m = 1, 2, 3, 4 modes as 
a function of time. Early in the simulation, during the for- 
mation of the ring, only the m = 1, 2 modes are important. 
Then all modes grow sharply when the ring breaks after 
t ~ 3500^ J , but m = 1,2 remain the most important ones, 
as the two-armed spiral seen in Fig 1 at t — 40017Q -1 suggests. 
After a couple of hundred binary dynamical times, however, 
the disc develops even finer spiral structure, which produces 
comparable amplitude for all Fourier modes in Fig. 5. 

We also measure the eccentricity of the disc, defined as 



\J d(f>i:(R,4>)v R e X p(i<l>)\ 



(2) 



Figure 6 shows the eccentricity profile both at the begin- 
ning of the simulation, and when we stopped it, 1200 time 
units later. Both profiles are built averaging over 20 snap- 



Figure 5. Evolution of the Fourier components of the disc Cm 
as a function of time. No mode dominates the structure of the 
disc after the initial transient period. 




Figure 6. Eccentricity of the disc as a function of radius, both 
at the beginning and the end of our simulation. The disc does not 
develop any significant eccentricity besides that due to the non- 
Keplerianity of the potential close to the binary, which is already 
present at the initial conditions. 



shots, to minimise the random fluctuations. The initial con- 
ditions of the simulation had the gas rotating in circular or- 
bits around the centre of mass of the system. However, the 
non-Keplerianity of the potential makes the gas rapidly ac- 
quire a small non-zero eccentricity, an effect that is obviously 
larger closer to the binary. After the collapse of the ring, the 
eccentricity of the disc reaches values e w 0.1, but then it 
decays again to smaller values. As shown in the figure, the 
eccentricity profile remains approximately the same as ini- 
tially until the end of the run, t « 1200f2 ( ^" 1 , even though the 
binary reaches e ~ 0.1 (Sect. 4). In a different run where the 
binary reached e ~ 0.3 the disc did not acquire a significant 
eccentricity either. These results are consistent with those of 
MacFadyen & Milosavljevic (2008), whose disc only became 
eccentric after a few thousand binary orbital times. 

To illustrate the long-term evolution of the disc shape, 
Fig. 7 shows a column density map averaging 100 snapshots 



at t 



700ft„ 



in a frame corotating with the binary. The 
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Figure 7. Density map averaged over tflo = 650—750, well after 
the initial transient evolution, plotted in a frame co-rotating with 
the binary. The disc shows no persistent asymmetry. 



image shows that the disc develops no persistent asymmetry, 
except for the streams taking material from the inner part 
of the disc to the black holes. 



4 EVOLUTION OF THE ORBIT 

The binary starts in a circular orbit of radius ao. However, as 
the binary exerts a positive total torque on the gaseous disc, 
it transfers its angular momentum to it, making the binary 
shrink. Less obviously, the coupling between the disc and 
the binary results in the excitation of binary eccentricity, 
even (as in our runs) in the case where the gas disc itself 
remains approximately circular. 



4.1 Analytic expectations 

The rate of angular momentum transfer between a binary 
and a surrounding accretion disc can be calculated analyti- 
cally (Pringle, 1991). Most relevant for our purposes is the 
work of Syer & Clarke (1995) and Ivanov et al. (1999), who 
calculated the evolution of a binary on a circular orbit, un- 
der the effect of an circum-binary a disc. They neglected 
the self-gravity of the disc, and assumed that the mass of 
the secondary black hole to be much smaller than the pri- 
mary. Although our conditions violate these assumptions, 
their calculations can provide us with an order-of-magnitude 
estimate of the effect we expect to obtain. 

In this section we follow Ivanov et al. (1999, their section 
4.1) to estimate the time-scale for binary shrinking. We first 
need to express the angular momentum transport in the disc 



due to self-gravity as a viscosity which is a power-law on 
surface density and radius, 



i/ocS it . 



(3) 



The parameters A; B are determined by the cooling mecha- 
nism in the disc. In our case the cooling time is proportional 
to the dynamical time-scale of the disc, t coo \ = f3/Q (Sec- 
tion 2), which sets the viscosity parameter to a = 0.4//3, 
(e.g., Rice et al., 2005). That, together with the condition 
Q = Q « 1.5 (Fig. 4), gives A = 2,B = 9/2, independent 
of the values of f3 and Qo- 

The other important parameter is 



S = 



27ra 2 So(a) 
W 2 ' 



(4) 



which is roughly the mass ratio between the disc and the sec- 
ondary black hole. £o is defined as the surface density the 
disc would have in the absence of the secondary. Over the 
relatively short time scales of the simulations self-gravity re- 
sults in relatively modest (less than an order of magnitude) 
changes to the disc surface density. To a reasonable approx- 
imation we can then extrapolate the E(i£) profile we use 
as initial conditions to R — a and get £o(a) ~ O.OlMa" 2 . 
That, together with the mass of the secondary, M2 = 0.25M, 
gives Sri 1/4. 

With the parameters obtained above, we follow Ivanov 
et al. (1999) calculation for the evolution of the surface den- 
sity of a circumbinary disc, under the assumption that the 
torque of the binary only affects a narrow ring in the inner 
part of the disc. We obtain that the time-scale for the black 
holes to merge is, for the case of a self-gravitating disc, 

1 M 3/2 



V0.93S 0.81aQ 2 7r 2 G 1 /2S2(a )a^ /2 

Im7 1 M 2 



(5) 



M 19.3aQ 2 G i/2£5/2 (ao)a 7/2- 

Using the parameters of our model, we estimate that tab 
3 x lO 4 ^ 1 . 



4.2 Numerical results 

Figure 8 shows the evolution of the binary orbital elements 
from our simulations. The orbit of the binary shrinks ini- 
tially at a rate da/dt ~ — 10~ 4 aof2o. The orbital decay gets 
somewhat stalled during the accumulation of material in a 
ring (t « 350S1Q 1 , see Fig. 1), but it recovers after the ring 
breaks. Practically the same rate of change is found for a 
binary that is initially in an orbit with e = 0.1. That par- 



ticular simulation was continued until t ~ 3300f2 







from 

that run it is clear that in the long term the rate of decay 
slows down, reaching da/dt ~ — 2 x 10~ 5 aof2o- We attribute 
this to the fact that the gas disc continuously absorbs angu- 
lar momentum from the binary and in average moves farther 
away from it - as this process goes on, it is increasingly hard 
for the gas to influence the binary. The order of magnitude 
of the decay rate is in good agreement with the theoretical 
expectations from the previous section 3 . 



3 Due to the long computing times required to complete these 
simulations, we have not yet explored the parameter space of dif- 
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Figure 9. Eccentricity as a function of semi-major axis for the 
different simulations, as shown in Fig. 8. Notice that the hori- 
zontal axis is inverted, as a decreases with time. The rate de/da 
remains roughly constant while the evolution of both parameters 
slows down with time. The thick lines show the results of a conver- 
gence test conducted with 8 million particles, the light curves refer 
to longer duration runs using the standard resolution of 2 million 
particles. The small amplitude variability results from the non- 
Keplerian potential produced by the disc, which means that the 
instantaneous Keplerian orbital elements oscillate slightly over 
the course of each orbit. 



Figure 8. Evolution of the binary orbit elements a (top panel) 
and e (bottom) for the different simulations. Solid lines show sim- 
ulations with the binary in an initially circular orbit, while dashed 
lines show binaries with e = 0.1 initially. In both cases the ec- 
centricity grows, and the orbital decay proceeds at the same rate. 
Thicker lines show shorter simulations with higher resolution that 
agree remarkably with the standard runs. 

We also study the evolution of the binary eccentricity. 
When the orbit of the binary is initially circular, it remains 
so during the first part of the calculation, during which the 
disc is still mostly symmetric (Fig. 5). However, after the 
disc becomes more asymmetric, the eccentricity starts grow- 
ing at a rate de/dt ~ 1.5 x 10~ 4 f2o- In a different simulation, 
where the binary orbit is initially eccentric, the eccentricity 
growth of the disc starts earlier, at t w , but pro- 

ceeds at approximately the same rate. As in the case of the 
binary size, the eccentricity growth slows down with time. 
It is interesting however to study the relation between ec- 
centricity and binary separation. Figure 9 shows that de/da 
remains approximately constant in the long simulation. We 
then conclude that the slower eccentricity growth at later 
times is mainly produced by the dearth of material to ab- 
sorb the binary angular momentum. It is not due to the onset 
of any intrinsic damping mechanism, such as might occur 
if resonant damping grows to balance resonant excitation 
of eccentricity (Moorhead & Adams, 2008). If eccentricity 
growth saturates in our case, it does so only for e £ 0.35. 

The standard resolution of 2 million particles was 
adopted based on a comparison between the physical torques 

ferent mass ratios. However, preliminary results show that lower 
disc masses result in slower binary decays, as expected. 



due to self-gravity and numerical torques due to code arti- 
ficial viscosity. It is not obvious that a resolution that min- 
imizes this numerical effect is also adequate to reproduce 
eccentricity growth accurately. Accordingly, we tested the 
numerical convergence of these results by running shorter 
simulations with 8 million particles, shown with thick lines 
in Figs. 8 and 9. These runs give the same results as our 2 
million particle runs. We are then confident that our sim- 
ulations capture properly the interaction between the disc 
and the binary. 

4.3 Maximum rate of orbital decay 

Ultimately, the main question that we wish to address is 
whether gas discs absorb binary angular momentum effi- 
ciently enough to drive the black holes toward coalescence on 
a short time scale. Evidently this will depend on the amount 
of gas that is typically present on very small scales within 
galactic nuclei. This quantity is not well known, either the- 
oretically or observationally. We can, however, estimate the 
maximum rate of orbital decay that can be driven by a sur- 
rounding thin accretion disc. Such a limit can be derived 
because, although more gas results in more rapid orbital 
decay, the amount of gas in the disc cannot be increased ar- 
bitrarily. Eventually, a very massive disc will fragment into 
stars rather than remaining as a fluid structure, destroying 
the hydrodynamic interaction that leads to shrinkage 4 . 

To estimate the maximum rate of decay we relax the 
constant /3 cooling time assumption that we introduced for 

4 The newly formed stars, of course, could refill the loss cone and 
thereby help to shrink the binary via stellar dynamical processes. 
We defer a calculation of that process to future work. 
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numerical convenience. In a real disc the cooling rate will 
be set by the physical conditions of the disc. Some discs 
will cool faster and fragment into stars, failing to trans- 
port angular momentum as a disc. We can then compute 
the maximum decay rate by finding the disc solution that 
maximizes \da/dt\, subject to the constraint that the disc 
is everywhere stable against fragmentation. We emphasize 
that the resulting disc model - which is marginally stable 
against fragmentation at every radius - is deliberately fine 
tuned and is not what one would expect as the outcome of 
gas inflow toward the nuclei. It suffices to yield a limit, but 
that limit is unlikely to be exactly realized in nature. 

From the analyses of Syer & Clarke (1995) and Ivanov 
et al. (1999), it is clear that the orbital decay will be faster 
for discs that are more massive (tdisc oc E _5//2 , see eq. 5). 
The first step is therefore to calculate the maximum surface 
density, as a function of radius, that a disc can have without 
fragmenting. We follow Levin (2007), and note that the frag- 
mentation boundary in effect represents the maximum stress 
that can be sustained in a quasi-equilibrium self-gravitating 
disc, and can be expressed in terms a maximum allowed 
value of a (Rice et al., 2005). If in addition the disc opacity 
is a function of temperature only, then there is an analytic 
solution for a critically self-gravitating disc with Q = 1. 
Each value of the temperature corresponds to unique values 
of SI and E, and together these specify the critical solution. 
We obtain 



time-scales for mass ratio q = 3 



8<7SB,2m p 4 6 T 



= 2m p ttGS 2 
fc b 1 Q ' ' 



k(T)E 
E 5 fr 7 , 



(6) 
(7) 
(8) 



where the opacity k(T) is taken from the Semenov et al. 
(2003) opacity tables. The calculation is stopped when the 
temperature reaches ~ 1300 K, as the opacity drops sharply 
there, allowing in principle unreasonably large disc masses. 
Moreover, the critical value of a can be much larger in this 
regime (e.g., Johnson & Gammie, 2003), rendering this anal- 
ysis invalid. 

Once a maximum E(i?) is established, we can calculate 
the maximum rate of binary shrinking using eq. 5. Figures 10 
and 11 show in solid lines this time-scale for systems with 
binary masses M = 3 x 10 5 ,3 x 10 6 ,3 x 1O 7 M , and mass 
ratios 1,3,9, as a function of binary separation. At large sepa- 
rations, the mass needed to make the disc fragment is small, 
and therefore the removal of angular momentum from the 
binary is not efficient. At short separations, much more mass 
is needed to make the disc fragment, so the process can be 
much faster. 

The sharp break that appears in every curve is due 
to the transition from the optically thin regime at large 
separations to the optically thick one. At large distances 
and low temperatures, the opacity is dominated by ice 
grains and goes roughly as k = k- 1cc T 2 , with Ki co ~ 2 x 
10" 4 cm 2 g _1 K- 2 (e.g., Bell & Lin, 1994). In addition, for 
a self-gravitating disc the disc properties are linked by 
T oc (E/S7) 2 , which gives an optical depth r oc E 5 /fl 4 . Look- 
ing at eq. 8, it is clear that the relation between surface den- 
sity and orbital frequency for a fixed a goes as E oc fi 11 / 10 for 
optically thin discs, while there is no constraint for optically 
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Figure 10. Time-scale for orbital decay for binaries with mass 
ratio q = 3 and different total masses. The solid line shows the 
maximum effect of a gaseous disc as found in this paper, while the 
dotted line shows the effect of stellar scattering from Milosavljevic 
& Merritt (2003) 



time — scales for M = 3e6 Msun 




distance [pc] 



Figure 11. As Fig. 10, but for binaries with total mass M 
3 X 10 Mq and different mass ratios. 



thick discs, as the E dependence cancels out. This explains 
the sharp decline of the time-scale at the point where the 
disc becomes optically thick - arbitrarily large masses are 
in principle allowed, making the angular momentum trans- 
fer through the disc very efficient. This behaviour, however, 
breaks down once the gas reaches T w 166 K and the opacity 
law changes to a shallower power-law on temperature. Fur- 
ther changes of slope are seen in the curves at even shorter 
separations, which are due to substructure in k(T). 

Overall, the time-scale for binary merger becomes 
shorter as the distance between the black holes is smaller, 
reaching a value of only ~ 3 x 10 8 yr for a 3 x 10 6 Mq binary 
at a separation of ~ 0.02 pc. The calculation is stopped at 
this point, as explained above, because the unperturbed disc 
would reach T » 1300 K. 

For comparison we also show in Figs. 10 and 11 the 
time-scale for orbital decay produced by scattering stars, 
tateiiar, as estimated by Milosavljevic & Merritt (2003), who 
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Figure 12. Top: 'Hung-up' time-scale for the evolution of the 
system including both stellar scattering and the effect of the disc. 
Bottom: Binary separation at that point. Binaries with masses 
M 10 7 Mq merge in less than a Hubble time. 



give 



isteiiar ~ 3 x 10 yr 



a 10 6 M 



© 



0.1a hard M 



(9) 



where ahard is the binary separation at the point where its 
orbital velocity is comparable to the galaxy velocity disper- 
sion a. For a, we take the value given by the M—a corre- 
lation, a « 60kms _1 (M/10 6 M ) ' 23 (e.g., Ferrarese et al., 
2006), which yields 



tstellar ~ 5 X 10 JT 



M 



(l + g) 2V 3x 10 6 M G 



>1.54 IPC 

a 



(10) 



Stellar scattering becomes very inefficient at sub-pc dis- 
tances from the binary, making the time-scale for decay 
longer than the Hubble time. The effect of the disc, in- 
stead, becomes more efficient at shorter distances. Where 
both time-scales intersect the evolution of the binary is the 
slowest. It is important to characterise this 'hang-up' point, 
as the time-scale there will give us an estimate of how fast 
the binary can merge, and whether it can do it in a Hubble 
time. Moreover, we would expect binaries to spend most of 
their evolution at that distance. Figure 12 shows the values 
of the hung-up time-scale t hu and separation a hu as a func- 
tion of binary mass M, for mass ratio q — 3 (there is no big 
difference in the range q — 1-10, as is clear from Fig. 11). For 
masses M 10 7 Mq , t nu is shorter than the Hubble time, 
therefore it is possible for such a binary to merge through 
the method described in this paper. Such binaries could most 
likely be caught at distances a nu ~0. 01-0. 06 pc. 

For convenience, we also present an analytical estimate 
for both ihu and dh u . From Figs. 10 and 11 it is clear that, for 
the masses we are interested in, the time-scales for both disc- 
driven and star-driven evolution intersect roughly where the 
disc-driven evolution turns due to the optically thin-thick 
transition. Numerically, we find that this occurs at r « 10. 
Using k oc T 2 , as appropriate for this regime, and eqs. 6-8, 
we find that the condition r = To translates to a fixed value 



4rom p 7rG, 



2/3 _ 



4<7SE 



1/3 



(11) 



= 7 .89xl0- 1 V 1 (^) 2 / 3 (^ ) - 8 / 3 (^) 1 /3 ; (i 2 ) 

which corresponds to a period of approx. 250 yr. The most 
common type of binary would then exhibit phase- dependent 
variations of the accretion rate on a time scale that is too 
long to be readily observable. Rarer binaries closer to merger 
would make better candidates for detection via accretion 
variations. A binary of total mass M will reach the hang-up 
point at a separation ahu = (GAf/li^) 1 ' 3 , 



M 1/3 
a hu(M) = — — -( 



< s.4/9^|16/9/ 9Q^ice ,2/9 

G 1 / 9 v 4r 7rm p ' Wo 1 4<j S b ' 



- 2 77 x 10~ 2 r 1 cf M ) 1/3 ( T ° i~ 4/9 C Qo ) w/9 ( - i 2/9 
-2.77X10 Pc( 3x106Mq ) ( 1Q ) ( 15 ) ( Q1 ) • 

We finally evaluate the stellar scattering time-scale at that 
position to obtain the hung-up time-scale, 

thu(Af) = 1.74 x 10 10 yr x (15) 

1 i M sl.21/70 n4/9/ Qo n-16/9/ a \-2/9 ,-, 

(l + g 2 ) l 3 x 10 6 M Q ' 4o j l iV l o.r ' 1 ' 

These analytical approximations are plotted as dashed lines 
in Fig. 12, and provide a good order-of-magnitude estimate. 



5 ACCRETION ONTO THE BLACK HOLES 

Gas that gets within O.lao of either black hole is taken away 
from the simulation to prevent the very short time-steps 
that would be necessary to follow it. This numerical trick is 
also useful to study the accretion onto each black hole, in 
particular its variability and behaviour as a function of the 
orbital phase. Figure 13 shows the accretion rate onto each 
black hole as a function of time. Overall, more gas is accreted 
by the secondary black hole, which is closer to the gaseous 
disc and has a greater specific angular momentum. Both 
curves show a large peak at the beginning, which is a result 
of the initial conditions - gas that was situated inside the 
tidal truncation radius. The accretion rates increase again 
at t ~ 400n,y 1 , after the disc loses its ring- like shape and a 
lot of radial mixing occurs. 

The accretion rates are strongly variable, especially 
on the time-scale of the orbital frequency. Once a M ~ 
3 x 10 6 Mq binary reaches a separation of 0.01 pc, the orbital 
frequency gets to an interesting observable range ~ 0.1 yr -1 . 
Observation of variability on these frequencies will thus help 
identify the presence of binaries in such close orbits, which 
would otherwise be extremely difficult to detect. One should 
note, however, that the maximal accretion rate that a self- 
gravitating disk can support is typically a small fraction of 



the Eddington limit '' . Additional gas infall to small radii 

5 The mass accretion rate for a critically fragmenting disc is M ~ 
37r 3 «QgG 2 (E(Q) /Q) 3 . In the optically thin regime this is only 
M ~ 10 — 6 Mq yr - 1 , roughly independent of the binary mass 
and separation. For smaller separations the accretion rate can 
grow larger as the disc becomes optically thick, but remains sub- 
Eddington for the whole range of our study. 
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Figure 13. Accretion rate onto each black hole as a function of 
time. 



(King & Pringle, 2007) - bypassing the region most vulner- 
able to star formation - would be needed to produce very 
luminous sources 

It should be noticed that the accretion rates are some- 
what influenced by numerical resolution. With 2 million SPH 
particles, an accretion rate of 10 -5 Mf2o corresponds to only 
100 accreted particles per orbital dynamical time. A shorter 
simulation with 8 million particles indeed shows up to factor 
2 discrepancies, however, the long-term trends and short- 
term variability are mostly the same. 



6 DISCUSSION 

In this paper we studied the interaction of a super- 
massive binary black hole with a self-gravitating circumbi- 
nary gaseous disc. Our focus has been on sub-pc scales, 
where binary hardening due to stellar dynamical processes 
is thought to be slow, and where efficient radiative cool- 
ing makes it likely that the gas will form a geometrically 
thin accretion structure. Studies of galaxy mergers, and in- 
direct evidence from our own Galactic Centre, suggest that 
gas flows to such small scales may be a common feature 
of galactic nuclei, particularly following mergers. If present, 
the existence of a gas disc around the binary can be a dom- 
inant factor in determining whether the binary merges and 
the time-scale for this process. The presence of a disc will 
influence the binary orbit, as it will transport the angular 
momentum of the binary outward. At the radii of greatest 
interest the evolution of such gas discs will be dominated 
by their self-gravity, which will alter the structure of the 
disc. Our simulations resolve angular momentum transport 
due to self-gravity physically, without requiring any ad hoc 
prescription to mimic viscosity. 

Our simulations show that the disc, after a transient 
phase due to the initial conditions, develops a quasi-steady 
spiral pattern that lasts for at least hundreds of dynamical 
times. The binary couples gravitationally to the disc, and 
loses angular momentum. This process makes the orbit of 
the binary decay at a rate da/dt ~ —few x 10~ 5 aofio (for our 
simulation parameters) , in good agreement with theoretical 



calculations for non-self-gravitating discs. The eccentricity 
of the binary also increases, at a rate de/dt ~ 10~ 4 fi . 

Due to the computational expense of the simulations, 
we have not followed the system for long enough to see a 
large change in the orbital parameters. Our direct numeri- 
cal experiments show only a modest (20%) decrease in the 
binary separation, accompanied by a growth of eccentric- 
ity to values of around e « 0.3. The evolution of both a, e 
appears to slow down with time, approaching asymptotic 
values a « 0.8ao;e w 0.4. However, we attribute this trend 
to the fact that the gaseous disc absorbs angular momen- 
tum and moves away from the binary, making the interac- 
tion less efficient. This is consistent with a viscous time-scale 
tviec ~ (R/H) 2 t OT b/a ~ SOOOfiy 1 at Ji ~ 3ao, where most of 
the disc mass is located early in the simulation. On the other 
hand, the ratio de/da remains roughly constant throughout 
our calculations, suggesting that the mechanism for eccen- 
tricity growth has not saturated. We speculate then that, 
given enough influx of external material (which could easily 
occur on the viscous time-scale described above), the orbit 
will continue its decay and the eccentricity will keep growing, 
perhaps reaching the very large values needed to influence 
the decay in the relativistic regime. 

We also studied the accretion onto the black holes dur- 
ing the orbital decay, and found that it is highly variable, 
especially on the time-scale of the orbital frequency. This 
variability pattern, if observed in AGN, could help us to 
identify binaries at sub-parsec separations. 

Motivated by the good agreement between the decay 
rate obtained from our simulations, and that predicted ana- 
lytically by Ivanov et al. (1999), we combine their formalism 
with that of Levin (2007) to calculate the maximum decay 
rate that can be obtained due to a self-gravitating disc. We 
find that the decay produced by the disc can dominate over 
stellar scattering once the binary separation is 0.01-0.1 pc. 
The time-scale for decay at the critical separation where gas 
disc and stellar processes are equally important is shorter 
than the Hubble time for binary masses M £ 10 7 Mq . This 
implies that geometrically thin gas discs could provide a so- 
lution to the 'last-parsec' problem for binaries in this range, 
whereas gas discs are unable to drive mergers within galaxies 
hosting the most massive black holes with masses of 10 s M Q 
and above. However, even for masses where gas discs can 
in principle result in mergers, the time scales are typically 
rather long. At low redshift, where the typical time scale 
between major mergers for galaxies hosting 10 Mq black 
holes is several Gyr (Volonteri et al., 2003), we would then 
expect that many galaxies would host binaries with sepa- 
rations close to the hang-up radius of a few hundredths of 
a parsec. At higher redshift the interval between mergers is 
much shorter - less than a Gyr - and interactions between 
three or more black holes in a single system could not be 
neglected. 

If we accept that most galactic mergers result in black 
hole mergers, the results presented here suggest that geo- 
metrically thin gas discs could in principle drive mergers on 
a short enough time scale for relatively low mass black holes 
at low redshift. Additional processes, either stellar dynam- 
ical or hydrodynamic, are needed if very massive binaries 
are to merge, and to avert ubiquitous 3-body interactions 
(that might result in ejections) at high redshift. One logi- 
cal alternative that we did not explore here is that provided 
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by even more massive discs. Such discs will fragment and 
form stars very close to the binary, and if those stars can 
be scattered by the binary prior to exploding as supernovae 
they will absorb the binary angular momentum almost as 
efficiently as if the disc had remained fluid. Given slower 
scattering, only low mass stars and stellar remnants will be 
able to contribute to binary decay. The stellar mass function 
is evidently critical in determining whether such a scenario 
is viable, as is the stellar dynamics of stars formed in a disc 
close to a (probably eccentric) binary black hole. 
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